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ABSTRACT 



We use deep far-infrared data from the FEP/GOODS-Herschel surveys and rest frame ultraviolet photometry to study the evolution of 
the molecular gas mass function of normal star forming galaxies. Computing the molecular gas mass, M„,oi, by scaling star formation 
rates (SFR) through depletion timescales, or combining IR luminosity and obscuration properties as in Nordon et al., we obtain My„„] 
for roughly 700, z = 0.2 - 3.0 galaxies near the star forming "main sequence". The number density of galaxies follows a Schechter 
function of M,noi. The characteristic mass M* is found to strongly evolve up to z ~ 1, and then to flatten at earlier epochs, resembling 
the infrared luminosity evolution of similar objects. At z~l, our result is supported by an estimate based on the stellar mass function 
of star foiTning galaxies and gas fraction scalings from the PHIBSS survey. We compare our measurements to results from current 
models, finding better agreement with those that are treating star formation laws directly rather than in post-processing. Integrating 
the mass function, we study the evolution of the M„,oi density and its density parameter Clmoi- 
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1. Introduction 

Stars form in cold, dense molecular clouds, and the molecular 
gas content of galaxies is an important constraint to galaxy evo- 
lution. In the interplay between accretion of gas, star formation, 
metal enrichment, and outflows, the molecular gas content re- 
flects the prevailing physical processes (e.g. Bouche et al. 2010; 
[Pave et a l. 2010; Lilly et al. 2013). For the cosmic baryon bud- 
get (e.g. Fukugita et al...l99 8), molecular gas can be of increas- 
ing relative importance at high redshifts where galaxies are gas- 
rich. 

Molecular gas is generally quantified using the CO molecule 
as a tracer. Local sam ples include several hundred targets (e.g. 
ISaintonge et alJ 1201 iL and references therein), but CO detec- 
tions of star-forming, intermediate/high redshift galaxies are still 
limited t o modest statistics of mostly very luminous galaxies 
(ICarilli & Walter 2013. and references therein). For normal star 
forming galaxies, the PHIBSS survey dTacconi et al.ll2013h de- 
rives scaling relations on the basis of CO detections of 52 nor- 
mal star forming galaxies at z ~ 1.2 and z ~ 2.2. In galaxies near 
the star-forming "main sequence" (M S, e.g. lNoeske et al.ir2007l; 
lElbazetal.ll2007t iDaddi et alJ l2007h . flie molecular gas mass, 
Mmoi, scales as SFR through a depletio n time scale Tiipp th at is 
only very slowly dependent on redshift (Tacco ni et al.ll2013l) . 



Keres et a l. (2003) reported the first attempt to derive the lo- 



cal molecular gas mass function, based on a sample of IRAS 
selected objects. To date, no derivation of a z > mass function 
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has been possible, due to the heterogeneity and paucity of the 
available detections. 

These difficulties have sparked interest in methods using dust 
mass as a tool to derive masses of the cold ISM of (distant) 
galaxies (e.g. Magdis et al. 20 llT 120121; FScoville 2012). Dust 
mass is converted into gas mass by adopting a metallicity, and a 
scaling of gas-to-dust ratio with metallicity. Typically, to derive 
dust masses, these methods require photometry on the restframe- 
submm tail of the SED, or accurate multi-band data closer to 
the rest far-infrared (FIR) SED peak. Such requirements are still 
not met individually for large samples of normal high-z galax- 
ies (Berta et al. in prep.). In a study of the restfr ame ultravio- 
let (UV ) and FIR properties of Herschel galaxies, iNordon et alJ 
(2013) established a method to derive M^oi of MS galaxies on 
the basis of their rest-FIR luminosity and rest-UV obscuration: 
properties which are readily available for larger samples. Here 
we apply the Tdep scalin g and the Nordon et al. recipe to the deep- 
est Herschel (Pilbratt et al]l2010h FIR extragalactic observations 
and related restframe UV data to derive Mmoi of MS galaxies, 
which are known to po wer ~90% of the cosmic star formation 
(iRodighiero et al.l201 ih . On this basis, we construct their molec- 
ular gas mass function and study its evolution since z ~ 3. 

Throughout this paper we adopt a ACDM cosmolog y with 
(/!,a„,nA) = (0.70.0.27.0.73) and a [Chabrieij (|2003[) initial 
mass function (IMF). 



2. Derivation of molecular gas mass and source 
selection 



We use the deepest PACS dPoghtsch et al.'20TO) far-i nfrared sur- 
vey, combining the PACS Evolutionary Probe (PEP, iLutz et all 
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1201 Ih and GOODS-Herschel (lElbaz et all r201lh data in the 
GOOp S-N and GOODS-S fields. We defer to 'Magne lli et alJ 
( 1201 3h for details about data reduction and catalog extra ction. 

An cilla ry data come from th e catalogs built by Be rta et all 
(1201 lb and lGrazianetal.1 (l2006i MUSIC), both including pho- 
tometry from the U t o Spitzer IRAC bands . To these we add 
Spitzer MIPS 24 mn ('MagnelH et al.l l201lh . IRS 16 fj.m data 
( lTeplitzetal.11201 1), GALEX D R6 data, and a coll ection of op- 
tical spectroscopic redshifts (see iBerta et all 1201 ll for details). 
When needed, photometric redshifts are used, reaching accura- 
cies of A(z)/(1 + z)=0.04 and 0.06 in the two fields, respectively 
(IBerta et al.il201 UlSantini et al.ll2 009). 

Based on current CO observations of local and z ~ 1 ob- 
jects (e.^. lSaintonge et aUlIoTTl 120121: pRcconi et aUlIoTsh . the 
relation between Mmd and SFR of MS galaxies can be described 
as a simple scaling with a depletion time scale mildly depen- 
dent on redshift, modulo uncertainties on the CO-to-molecular 
gas mass conversion factor aco'. M,noi/SFR = Tdep, with Tdep - 
1.5 10^x(l -Hz)"' [Gyr]. Second ordereffects, likely driven by the 
actual distrib ution of dust and gas in galaxies, have been stud- 
ied bv Nordon et al.l (1201 3h . while investigating the infrared ex- 
cess (IR X) and the observed UV slope /? of 1.0 < z < 2.5 PEP 
galaxies (iMeurer et al.lll999l) . Using a smaller sample of z > 1 
galaxies with CO gas masses, Nordon et al. study the link be- 
tween AiRx and molecular gas content, finding a tight relation 
between the scatter in the Airx-/3 plane and the specific attenua- 
tion contributed by the molecular gas mass per young star. Herqj 
AiRx = 2.51og(SFRiR/SFRuv + 1) represents the effective UV 
attenuation. These authors derive the equation 



On the FIR side, we apply a 3 <t flux cut at 160 //m, the band 
that best coiTelates with L(IR) ("Elbaz et al]l201ll: iNordon et aP 
2012). When deriving UV parameters, it is important to avoid 
contamination of observed bands by the 2100 A carbonaceous 
absorption feature. Combining GALEX and optical photometry, 
we define four redshift windows: 0.2 < z < 0.6, 0.7 < z < 1.0, 
1.0 < z < 2.0, and 2.0 < z < 3.0. This choice reduces the loss 
of sources due to restframe UV requirements to <2%. The total 
number of sources in each redshift bin is included in Table[T] Our 
sample contains 43 galaxies hosting an X-ray AGN, of which 
only 4 are Type- 1 . Distributions of L(IR) and /3 for the AGN 
galaxies are similar to those of inactive ones. We have verified 
that our conclusions below are not affected by the inclusion of 
these AGN hosts. 

3. The molecular gas mass function 

The comoving number density of galaxies in intervals of Mmd is 
computed adopting the well known l/Va formalism: 



<D(M)AM ^Yj—AM 



log 



Mrr 



lilO^Mo] 



= l°g| 7^^-^ 1-0.20 (AiRx - 1.26yS)+0.09,(l) 

[Moyr-i' 



where the SFR includes both IR and UV contributions. This pro- 
vides an estimate of molecular gas mass (already including the 
contribution of Helium) using only widely available integrated 
rest-frame UV and FIR photometry, consistent with z ~ 1 deple- 
tion times by calibration. We use the FIR and UV photometry to 
apply the Tacconi et al. (2013) rjep scaling and the Nordon et al. 
(1201 3h recipe to ~ 700 MS galaxies i n our survey. 

SFRiRuv are computed using the iKennicuta (Il998h calibra- 
tions. In frared lu minosities, L(IR), are derived by SED fitting 
with the Berta et al. (2013) templates hbrary. A lternative meth- 
ods (e.g. Wuvts e t al. 2011; Nor don et alJ2012j) lead to ea uiva- 
lent results within a 10-20% scatter Csee iBerta et al.ll2013l) . UV 
parameters are derived from the restframe 1600 and 2800 A lu- 
minosities, following the procedure of Nordon et al. (2013). 

Statistical uncertainties on M^oi have been computed per- 
forming 10000 Monte Carlo realizations while varying photom- 
etry and redshifts according to their observed uncertainties and 
assuming Gaussian error distributions. Including also the intrin- 
sic scatter of Eq. [1] the median total uncertainty on M^oi is 
~40%, and only roughly 10% of the sample have mass uncer- 
tainties larger than 50%. 

Sample selection is driven by the need to compute total SFRs 
and by the requirements imposed by Eq. [1] First of all, we limit 
the analysis to galaxies lying within |Alog(SFR)Ms| ^ 0.5 from 
the MS. The MS is assumed to have unit slope in the stellar mass 
vs star formation rate (M*-SFR) plane, and a specific-SFR nor- 



malization varying with redshift as sSFRms [Gyr ] 
(lElbazetal.1120111) . 



26xr 



-2.2 



(2) 



where the sum is computed over all sources in the given Mmoi 
bin. The accessible volume V^ is a spherical shell delimited by 

z|™ and min(zmax,Zmax)- Here Zmax is the maximum r edshift at 
which a galaxy would be observable in our survey (ISchmidtl 
'1968'), and z^'" define redshift bins. 

-^' min.max 

Source selection is mainly based on a 160 fim cut, and the 
UV requirements do not produce significant source losses. On 
the other hand, the conversion from 5(160) to L(IR) is obtained 
adopting a family of SED templates spanning a variety of colors. 
Moreover UV properties are also involved in computing total 
SFR and in applying Eq. 1 . 

As a consequence, the conversion between 160 yum fluxes 
and molecular gas mass is not unique, but comprises a distri- 
bution of "mass to light ratios". Because of this scatter in Mmoi 
vs. 5(160), a flux cut induces a molecular mass incompleteness. 
This effect was thoroughly studied by Fontana et al. (2004) and 
Be rta et al.l (1^07), among others, in the case of stellar mass. 
Equivalently, the recipe defined by these authors can be applied 
to the specific case of molecular gas mass and FIR fluxes to de- 
rive completeness corrections as a function of Mmd. Using the 
distribution of the parent MS population, and comparing it to 
that of PACS galaxies leads to similar results. 

The comoving number density of galaxies is shown in Fig. 
[U Results from the Tdep scaling or Eq. [1] are very similar; thus 
only the latter are shown. Table [T] reports both estimates, along 
with average completeness values for each mass bin, which also 
includes the photometric completeness of the FIR catalogs. 

An independent characterization of the mass function is 
provided by the maximum likelihood approach by Sandage, 
Tammann & Yahil (1 19791 STY) . We adopt the Bayesian imple- 
mentation by Berta et al.l (l2007h and adapt it to our case, thus 
fully propagating the actual Mmoi uncertainties into the paramet- 
ric function evaluation. The adopted functional form to describe 
the mass function is a Schechter ( 1976.) function: 



0(M)£/M ■■ 









e "• dM, 



(3) 



' SFRiR is based on bolometric IR luminosity; SFRuv is derived from 
UV continuum, uncorrected for dust attenuation. 



where <!)* represent the normalization, a the slope in the low 
mass regime, and M* the transition mass between a power-law 
and the exponen tial drop-off and th e e-folding mass of the lat- 
ter. We defer to iBerta et aP (l2007h for more details about this 
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Fig. 1. The molecul ar gas mass functi on of Hersch el galaxies. Left: co mparison of t he l/Va estimate (green squares, bas ed on Eq. 
[B to Uterature data (iKeres et al.ll2003h and models JObreschkow & Raw lings 20 0931; iLagoset alj|201 ulPuffv et al.ll2012h . The red 



line and shaded area at z- 1.0-2 .0 are obtained scaling the Muzzi n et al.i (i2013i) z- 1 .0 - 1 .5 stellar mass function using molecular 
gas fractions derived bv lTacconi et al. (2013). When needed, masses found in the literature were scaled by the factor 1.36 necessary 
to account for Helium, and matched to our set of cosmological parameters. Right: results of the parametric STY evaluation of a 
Schechter mass function (red lines) and its 3 cr uncertainty (shaded areas). 



method. Table |2] includes the most probable parameter values 
and their 3 cr uncertainties, obtained with both M,™! estimates. 
The right-hand panel of Fig. [T] compares the result of the STY 
analysis to the 1 /Va mass function, based on Eq.[T] No STY anal- 
ysis is attempted for the highest redshift bin, because we cover 
only the very massive end at these epochs. 

The Schechter M* parameter is found to increase by more 
than a factor of 3 between z - 0.4 and 0.8, and then to 
flatten at z > 1 . This rate resembles the evolution found by 
iGruppioni et al.l (1201 3l) in the IR luminosity function of normal 
star forming galaxies, and reflects the link between SFR and 
Mmoi- At the same time, we find that <!)* varies by only a fac- 
tor of 2 over the 0.4-1.5 redshift range. The net elTect is a sig- 
nificant evolution of the number density of galaxies with large 
^moi, while at the low mass end it remains roughly constant. 
Finally, the most probable value of a steepens as redshift in- 
creases, but this might be simply an effect of the different mass 
ranges effectively constrained at different redshifts (note that the 
3 cr confidence levels are consistent with nearly no evolution). 

We compare this Mmoi mass function to an estimate based 
on ste llar mass. Using the PHIBSS CO survey, Tacco ni et al. 
(1201 3L see their Fig. 12) compute molecular gas fractions for 
Z = 1 .0 - 1 .5 normal star forming galaxies, as a function of M* . 
We combine these gas fractions with t he z = 1 .0- 1 .5 stella r mass 
function of star forming galaxies by Muz zin et al.l (1201 3h. who 
separate star forming an d quiescent systems (see also lllbert et alJ 
l2013l:lDrorv et al.ll2009h . To account for the MS widfli in this es- 
timate, we apply a 0.2 dex Gaussian smoothing. The result is 
close to the observed mass function at z = 1.0 - 2.0 (Fig.[T]i. 



4. Discussion 

Our observed Mmoi function is computed only for normal star 
forming galaxies within +0.5 dex in SFR from the main se- 



quence. Thus it represents a lower limit to the total M,noi mass 
function, by missing passive (low SFR) galaxies and powerful, 
above-sequence objects. Passive galaxies will provide a neg- 
ligible contribution unless they include a hypothetical popu- 
lation of molecular gas-rich galaxies forming stars with very 
low efficiency. For reference, local passive galaxies, lying at 
A(sSFR)ms ^ -1.0 dex , have molecular gas fractions /moi < 2% 
(ISaintonge et al.ll2012h . Also star bursting galaxies, i.e. those ly- 
ing above the MS , should play on l y a rn inor role in shaping 
the mass function. i Rodighiero et al.l (1201 ll) have shown that the 
contribution of above-sequence sources to the number density 
of star forming galaxies and total star formation rate density at 
z - 1.5 - 2.5 is small. For our A log(SFR)MS - 0.5 cut, we find 
that objects above the MS contribute no more than ~ 4% to the 
numberdensity of star fomiing galaxies, and ~ 16% to their SFR 
density. Depending on whether galaxies rise above the main se- 
quence mostly due to an increased star formation efficiency at 
fixed gas mass, or due to larger gas masses, the effect on the 
mass function will vary between a global upward ~ 4% shift 
or a preferential increase at higher gas masses within the limits 
permitted by the 16% SFR contribution. Saintonge et al. ( 2012) 
show that in local starbursts the two effects share a 50%-50% 
role in causing SFR changes with respect to the MS. 

Ou r results are compared in Fig. [T] to the iKeres et al] 
(l2003h local CO mass function of IRAS-selected galaxies. 
Applying a variable CO-to-H2 conversion factor, inversely de- 
pendent on L(CO) itself, Obreschkow & Rawlings ( 2009b a) ob- 
tained a revisited mass function quickly dropping at the high- 
mass end. Modehng was implemented by post-processing the 
iDe Lucia & BlaizoB(l2007h semi analytic (SAM) results, and as- 
signing the atomic and molecular gas content to galaxies via a 
set of physical prescriptions (see Obreschkow et al. 2009). Their 
expectations (Fig. [1]) tend to predict a too steep mass function 
both at low and high redshift. 
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^^^' . The molecular gas density parameter peaks at z ~ 
1.0; for comparison Ohi, derived from damped Ly-a systems 
(Prochaska & Herbert-Fort 2004; Rao et al. 2006), and includ- 
ing also low column density cases, does not evolve between 
z=0.5and4.0. 

Combining the deepest FIR Herschel extragalactic observa- 
tions available and restframe UV information for roughly 700 
main sequence galaxies, we have built the first molecular gas 
mass function at redshift z > 0, and the first estimate of its den- 
sity evolution up to z = 3. While future mm/submm surveys 
will significantly improve our knowledge of the evolution of the 
molecular content in galaxies, we have provided a basis for re- 
finement of galaxy evolution models that are accounting for the 
molecular phase. 
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Fig. 2. 7b;?: redshift evolution of the molecular gas mass den- 
sity, based on Eq. [1] Red symbols and solid error bars be- 
long to the total mass density; dashed error bars mark lower 
limits limited to the mass range covered by observations. The 
black star is comp uted by integrating the local mass function 
dKeres et al]|2003l). Grey squa r es and tria ngles represen t the to - 
tal M* density bv lllbert et all (l2013l) and lMuzzinetal.l(l2013h . 
respectively. Light blue symbols belong to star forming galaxies 
only. Bot tom: evolution of Qmoi (re d symbols, th is work) a nd Qhi 
(Prochas ka & Herbert-F ort 2004; IZwaan et all 2005; R ao et all 
12006 ). where we have divided the latter by 1.3 in order to avoid 
counting Helium twice. 



A sec ond model based o n the SAM approac h was devel- 
oped by iLagos et aP (1201 ih . starting from the iBower et al.l 
(2006) galaxy formati on mod el. As in the case of 
[Obresch kow & RawUngsl ( l2009ah the iBUtz & Rosolowskvl 
(.2006) star formation law is adopted, but in this case it is im- 
plemented throughout the galaxy evolution process in the SAM 
model. The model was tested against the observed stellar mass 
density evolution and the atomic and molecular gas content of 
local galaxies. Results are shown in Fig. [1] and are now much 
closer to our observed mass function up to z ~ 2. 

Figure [U fi nally includ e s pred ictions of the hydrodynamical 
simulation by iDuffV et al] (I2012D . which overall tend to over- 
estimate the mass function, but are consistent with data at the 
observed high-mass tail. 

Adopting the STY Schechter results, we integrate the mass 
function from 10^ M© to infinity and obtain a measure of the 
molecular gas mass density (see Table [2] and top panel of Fig. 
|2]l. A second estimate obtained by limiting the integral to the 
mass range effectively covered by observations is also provided. 
The molecular gas mass density increases by a factor of ~ 4 
from z = to z = 1 and then remains almost co nstant up to 
z ~ 2 . For reference we plot also the M* density ( Muzzinet alJ 
1201 3t lllbert et all l20T3h up to z = 4. The diff^erent trends re- 
flect the growth of gas fractions as a function of redshift (e.g., 
lTaccometaT]l2013 !). 

Finally, we compute the redshift evolution of the density 
parameter i^moKz) - Pmo\{z) I pdz) (Fig. |2] and Table lU, where 
the critical density at the given redshift is given by pdz) - 
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Table 1. The molecular gas mass function of Herschel galaxies, derived with the 1 /Va method. 







0.2 < c < 0.6 






0.7 <z< 1.0 






1.0<z<2.0 
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Eq.Q 
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Eq.Q 
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I"dL-p 




[A7I Mo] 


[10-*ft^„Mpc-'dex-'] 




[10-^ /,^„ 
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9.0 


23.42±6.26 


- 


0.52 


- 


- 


- 




- 


- 


- 


- 
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9.2 


25.38±4.80 


18.42±5.11 


0.67 


- 


- 


- 


- 


- 


- 


- 


- 


- 


9.4 


22.01 ±4.02 


24.42±4.39 


0.78 


- 


- 


- 


- 


- 


- 


- 


- 


- 


9.6 


13.00±2.77 


17.86±3.32 


0.86 


12.84±3.87 


20.15±6.08 


0.61 


- 


- 


- 


- 
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- 


- 
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Table 2. Results of STY analysis, and molecular gas mass density. Top: results based on Eq. [T] Bottom: results obtained with the 
Tdep scaling of SFR. 
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0.2 < 
Most piob. 


z < 0.6 
Min. 
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0.7 <z< 1.0 
Most prob. Min. 
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a 


11.0 
10.34 
^0.91 


9.8 

10.23 
-1.08 


11.0 
10.34 
-0.90 


7.8 
10.84 
-1.12 


7.4 
10.72 
-1.15 


10.6 
10.86 
-1.03 


5.6 
10.89 
-1.16 


5.3 
10.84 
-1.16 


6.6 
10.94 
-1.04 


- 


- 


- 


p„„i(tot) [10' A70 Mo Mpc"-'] 
p„,„,(lim)[10'/i7(,MoMpc-^] 


2.266 
2.166 


1.751 
1.619 


2.289 
2.188 


5.974 
4.173 


4.317 
3.231 


7.818 
5.361 


4.836 
3.958 


4.112 
3.301 


5.711 
4.972 


1.533 


1.184 


1.882 


n™,(tot)[10-'A7^] 
n„K,i(lim)[10-5/,-^] 


11.34 
10.84 


8.76 
8.10 


11.46 
10.95 


18.00 
12.57 


13.01 

9.74 


23.56 
16.15 


7.19 
5.89 


6.12 
4.91 


8.50 
7.40 


0.92 


0.71 


1.13 




Based on T^^p 


0.2 < 
Most piob. 


z < 0.6 
Min. 


Max. 


0.7 <c< 1.0 
Most prob. Min. 


Max. 


1.0 
Most prob. 


< z < 2.0 
Min. 


Max. 


2.0 < 
Most prob. 


z<3.0 
Min. 
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<I>* [10-*/!^ Mpc-'dex"'] 

log(M*) [h-;^ Mo] 

a 


12.0 
10.32 
-1.01 


9.8 

10.27 
-1.11 


13.6 
10.36 
-1.00 


7.4 
10.79 
-1.06 


6.8 
10.70 
-1.10 
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10.83 
-1.00 
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10.87 
-1.11 


4.3 
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-1.04 
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- 
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1.963 
1.810 
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3.437 
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4.407 
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2.590 
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12.01 
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9.06 
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10.97 
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13.28 


5.77 
4.81 


4.77 
3.85 


8.04 
6.99 


0.73 


0.56 


0.90 



Notes. Minimum and maximum values are computed at the 3 cr confidence level. The molecular gas mass density p(tot) is obtained by integrating 
the mass function between 10' M© and infinity; p(lim) is computed solely on the mass range covered by PEP data (see Table[TJ. 






